# -*- coding: utf-8 -*-
"""
Created on Tue Oct 01 10:12:42 2013
field performance
@author: Maxim
"""
import aircraft
from performance import ExtendedInput, flightMechanics
from numpy import linspace, zeros, sin,cos,tan, arcsin, ones,cumsum

class _FieldRslt:
    """
    Stores and displays result of field simulation
    """
    def __init__(self,name=None,n=3):
        self.name = name
        self.fuelBurned = zeros(n)
        self.distance = zeros(n)
        self.time = zeros(n)
        self._i = 0

    def display(self, detailed=False):
        if not self.name==None:
            print '==>',self.name
        print 'Fuel burned = ',self.fuelBurned
        print 'Time = ',self.time
        print 'Distance = ',self.distance
        
    def append(self,distance,time,fuel):
        i = self._i
        self._i += 1
        self.distance[i] = distance
        self.time[i] = time
        self.fuelBurned[i] = fuel
    
    def _summarize(self):
        self.time = sum(self.time)
        self.distance = sum(self.distance)
        self.fuelBurned = sum(self.fuelBurned)

class Field:
    def __init__(self,aircraft,fieldAltitude=0, hObstacle=15):
        self.ac = aircraft
        self.tm = aircraft.analysis.thrust
        self.hField = float(fieldAltitude)
        self.hOb = float(hObstacle)
        self.fm = flightMechanics()
        self.bi = ExtendedInput(self.ac,altitude=self.hField)
    def takeoffSimulation(self,powerSetting=95.0,CLto=0.7,CLmaxTo=1.9,muClean=0.02,
                   dtRotation=1.5):
        """
        Takeoff simulation. Returns data structure containing information about 
        takeoff path.
        
        Parameters
        ----------
        
        CLto : float
            lift coefficient at ground roll
        CLmaxTo : float
            maximum lift coefficient in takeoff configuration
        muClean : float
            landing gear friction coefficient in clean configuration (no brakes)
        """
        self._get_missing_data(CLmaxTo)
        rslt = _FieldRslt('Take off',3)
        Sgr, tgr, wfgr = self._TO_groundroll(CLto, powerSetting, muClean)
        rslt.append(Sgr, tgr, wfgr)
        Srot, trot, wfrot = self._TO_rotation(dtRotation)
        rslt.append(Srot, trot, wfrot)
        Sco, tco, wfco = self._TO_climbout(wfgr+wfrot, CLmaxTo, powerSetting)
        rslt.append(Sco, tco, wfco)
        return rslt
    
    def landingSimulation(self,CLgr=0.7,CLmaxLDG=2.183,muBreaks=0.3,idlePower=0.0,
                   loadFactor=1.2,reactionTime=1.0):
        """
        Landing simulation. Returns data structure containing information about 
        landing path.
        
        Parameters
        ----------
        
        CLgr : float
            lift coefficient at ground roll
        CLmaxLDG : float
            maximum lift coefficient in landing configuration
        muClean : float
            landing gear friction coefficient with brakes
        idlePower : float, Watt
            idle power at approach in watts. 0 value is preferred due to 
            temporary analysis limitations
        loadFactor : float
            load factor at approach speed
        """
        self._get_missing_data(CLmaxLDG)
        rslt = _FieldRslt('Landing',2)
        Sap, tap, wfap = self._LDG_approach(idlePower,0.0,loadFactor)
        rslt.append(Sap,tap,wfap)
        Sgr, tgr, wfgr = self._LDG_groundroll(0.0, CLgr,muBreaks,reactionTime)
        rslt.append(Sgr, tgr, wfgr)
        return rslt
    
    def _get_missing_data(self,CLmax):
        # ground effect calculation
        hw = self.ac.wing.aapex[2] - self.ac.landingGear.groundContact_Z[2]
        b = self.ac.wing.span
        _hw16 = (16.0*hw/b)**2.0
        self.Ge = _hw16 / (1.0+_hw16)
        #speeds
        self.bi.set_CLmax(CLmax)
        self.Vstall = self.fm.vStall(self.bi)
        self.Vlo = 1.1 * self.Vstall
        self.Vtr = 1.2 * self.Vstall
        self.Vflare = (self.Vlo + self.Vtr) / 2.0
        self.Vap = 1.3 * self.Vstall
        self.Vtd = 1.1 * self.Vstall
        self.Va = (self.Vap + self.Vtd) / 2.0
        # weight
        self.W = self.ac.mass() * self.bi.g
    
    def _TO_groundroll(self, CLto, powerSetting=95.0, mu=0.02, nSeg=10):
        tol = 1e-5
        err = tol+1.0
        T  = zeros(nSeg)
        ff = zeros(nSeg)
        dt = zeros(nSeg)
        V = linspace(0.0,self.Vlo,nSeg)
        CDto = self.bi.Cd0 + self.Ge*self.bi.k*CLto**2.0
        dV = V[1]-V[0]
        Q = 0.5 * self.bi.density * V**2.0
        L = Q*self.bi.refArea*CLto
        D = Q*self.bi.refArea*CDto
        for i,v in enumerate(V):
            self.tm.runAnalysis(v,self.bi.density,powerSetting)
            T[i] = self.tm.thrust
            ff[i] = self.tm.fuelFlow
        #FIXME: thrust at V=0 fails to converge due to propeller analysis limitation
        T[0] = T[1]
        while err>tol:
            tLast = sum(dt)
            Wf = self.bi.g*dV*cumsum(ff*dt)
            W = self.bi.wt*ones(nSeg)-Wf
            F = T-D-mu*(W-L)
            dt = (W/self.bi.g)/F
            err = abs(tLast-sum(dt))
        ds = V*dt
        t = cumsum(dt*dV)
        s = cumsum(ds*dV)
        self.ffVlo = ff[-1]
        return s[-1], t[-1], W[0]-W[-1]

    def _TO_rotation(self,dtRot=1.5):
        Srot = self.Vlo*dtRot
        wf = self.bi.g * self.ffVlo*dtRot
        return Srot, dtRot, wf

    def _TO_climbout(self, fuelBurned, CLmax, powerSetting=95.0):
        """ climbout path of takeoff
        """
        Wtr = self.bi.wt - fuelBurned
        nFlare = self.bi.density*self.Vflare**2/2.0 * self.bi.refArea*CLmax / Wtr
        R = self.Vflare / nFlare
        self.tm.runAnalysis(self.Vtr,self.bi.density,powerSetting)
        fuelFlow = self.tm.fuelFlow
        Q = self.bi.density*self.Vtr*self.Vtr / 2.0
        CL = Wtr / (Q*self.bi.refArea)
        CD = self.bi.get_drag_coefficient(CL)
        D = Q*self.bi.refArea*CD
        climbAngle = arcsin((self.tm.thrust-D)/Wtr)
        #-->
        Sflare = R*climbAngle
        hFlare = R*(1.0-cos(climbAngle))
        tFlare = Sflare / self.Vflare
        WfuelFlare = fuelFlow*tFlare
        #-->
        hTR = self.hOb - hFlare
        RCtr = self.Vtr*sin(climbAngle)
        ttr = hTR/RCtr
        Str = self.Vtr*cos(climbAngle)*ttr
        WfTR = self.bi.g*fuelFlow*ttr

        Sclimb = Sflare + Str
        timeClimb = tFlare + ttr
        WfClimb = WfuelFlare + WfTR
        return Sclimb, timeClimb, WfClimb
        
    def _LDG_approach(self,idlePower,idleFuelFlow, loadFactor):
        T = idlePower * self.Va / 2.0
        ff = idleFuelFlow
        CL = self._get_CL(self.Va)
        CD = self.bi.get_drag_coefficient(CL)
        theta = arcsin(1.0/(CL/CD)-T/self.W)
        R = self.Va**2.0 / (self.bi.g*(loadFactor - 1.0))
        hf = R*(1.0 - cos(theta))
        Sa = (self.hOb - hf) / tan(theta) + R*sin(theta)
        Vg = self.Va * cos(theta)
        ta = Sa / Vg
        Mfa = ff * ta
        return Sa, ta, Mfa*self.bi.g

    def _LDG_groundroll(self,idleFuelFlow,CLgr,mu,reactionTime,nSeg=10):
        CDldg = self.bi.Cd0 + self.Ge*self.bi.k * CLgr**2.0
        dV = -self.Vtd / (nSeg-1)
        V = linspace(self.Vtd,0.0,nSeg)
        Q = self.bi.atm.density*V*V / 2.0
        L = Q*self.bi.refArea * CLgr
        D = Q*self.bi.refArea * CDldg
        tol = 1.0e-6
        err = tol+1.0
        dt = ones(nSeg)
        ff = ones(nSeg)*idleFuelFlow
        while err>tol:
            tLast = sum(dt)
            Wf = self.bi.g * cumsum(ff*dt)
            W = self.bi.wt*ones(nSeg) - Wf
            F = -D - mu*(W-L) #ma
            a = F * self.bi.g/ W
            dt = dV / a
            err = abs(tLast-sum(dt))
        ds = V*dt
        t = cumsum(dt)
        s = cumsum(ds)
        Sgr = s[-1] + reactionTime*self.Vtd
        tgr = t[-1] + reactionTime
        Wfgr = Wf[-1]
        Mfgr = Wfgr / self.bi.g
        return Sgr,tgr,Mfgr*self.bi.g
    
    def _get_CL(self,V):
        return 2.0*self.W / (self.bi.density*V*V*self.ac.wing.area)


def run_field(aircraft,fieldAltitude,obstacleHeight=15.,takeoffPowerSetting=100.,
              CLtakeoff=0.7,CLmaxTakeoff=1.9,muClean=0.02,reactionTimeTakeoff=1.5,
              CLgroundroll=0.7,CLmaxLanding=2.2, muBreaks=0.3,idlePower=0.0,
              landingLoadFactor=1.2,reactionTimeLanding=1.0):
    """
    Runs takeoff and landing simulation.
    
    Parameters
    ----------
    aircraft : object
        aircraft object
    fieldAltitude : float, m
        altitude of takeoff and landing field
    obstacleHeight : float, m
        height of the obstacle
    takeoffPowerSetting : float
        engine power setting at takeoff [0,100]. High power setting is preferred.
    CLtakeoff : float
        takeoff ground roll lift coefficient
    CLmaxTakeoff : float
        maximum lift coefficient of aircraft in takeoff configuration
    muClean : float
        landing gear friction coefficient in clean configuration
    reactionTimeTakeoff : float, sec
        pilot reaction time at takeoff
    CLgroundroll : float
        landing ground roll lift coefficient
    CLmaxLanding : float
        maximum lift coefficient of aircraft in landing configuration
    muBreaks : float
        landing gear friction coefficient with breaks
    idlePower : float
        power setting when engine is idle
    landingLoadFactor : float
        aircraft load factor at landing
    reactionTimeLanding : float, sec
        pilot reaction time at landing
    """
    field = Field(aircraft,fieldAltitude,obstacleHeight)
    takeoff = field.takeoffSimulation(takeoffPowerSetting,CLtakeoff,CLmaxTakeoff,
                                      muClean,reactionTimeTakeoff)
    landing = field.landingSimulation(CLgroundroll,CLmaxLanding,muBreaks,idlePower,
                                      landingLoadFactor,reactionTimeLanding)
    takeoff._summarize()
    landing._summarize()
    return takeoff, landing


def run_test1():
    ac = aircraft.load('V0510')
    to, ldg = run_field(ac,0,15)
    to.display()
    ldg.display()

if __name__=="__main__":
    run_test1()